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We propose two methods for generating non-Gaussian maps with fixed 

power spectrum and bispectrum. The first makes use of a recently proposed 
rigorous, non-perturbative, Bayesian framework for generating non-Gaussian 
distributions. The second uses a simple superposition of Gaussian distribu- 
tions. The former is best suited for generating mildly non-Gaussian maps, and 
we discuss in detail the limitations of this method. The latter is better suited 
for the opposite situation, i.e. generating strongly non-Gaussian maps. The 
ensembles produced are isotropic and the power spectrum can be jointly fixed; 
however we cannot set to zero all other higher order cumulants (an unavoid- 
able mathematical obstruction). We briefly quantify the leakage into higher 
order moments present in our method. We finally present an implementation 
of our code within the HEALPIX package Q. 



I. INTRODUCTION 

In most inflationary scenarios the fluctuations in the matter fields generated by the 
oscillating inflaton display Gaussian statistics. This requires the temperature fluctuations in 
the Cosmic Microwave Background (CMB) to be Gaussian distributed to a very high degree 
of accuracy at sufficiently large angular scales. On smaller scales the effect of late time non- 
linear evolution will introduce a certain amount of non-Gaussianity in the CMB. The study 



of the non-Gaussianity of CMB fluctuations is therefore crucial to both the understanding 
of the fundamental processes generating the fluctuations and to the understanding of the 
various foreground and astrophysical contributions. 

In addressing an issue of this nature with reference to observational strategies it is im- 
portant to be able to simulate CMB maps with non-Gaussian signatures. These can then be 
used in the refinement of estimation techniques and the design of evermore accurate satellite, 
balloon-borne, and ground-based experiments. It seems therefore desirable to develop fast 
algorithms for simulating not only Gaussian signals (as extensively done in the past, e.g. 
HI]), but also maps allowing for non-Gaussianity. In the past the bispectrum has proved an 
invaluable tool in studying CMB non-Gaussianity (see for instance [0-0). Also algorithms 
generating Gaussian maps usually use the power spectrum as a controlling parameter. We 
therefore seek to complement earlier work by producing an algorithm for generating non- 
Gaussian maps with fixed power spectrum and bispectrum. As a matter of fact the algorithm 
we shall propose can be extended to fix simultaneously any higher order moment, naturally 
at a computational cost. 

The difficulty in generating non-Gaussian maps is in part due to the lack of suitable 
Probability Distribution Functions (PDFs) with the required parametrisation, and in part 
due to the requirement that non-zero higher moments require the multivariate generation 
of correlated sets of modes. The latter is imposed by statistical isotropy, in all cases except 
in Gaussian maps. We shall address the first of these problems with reference to the re- 
cent work in 0, in which a rigorous, non-perturbative, Bayesian framework for generating 
non-Gaussian distributions was proposed. We also propose an alternative PDF, a simple 
superposition of Gaussian distributions. We then address the second problem with a rather 
simple trick for creating the required mode correlations enforcing isotropy. 

This paper is organised as follows. In section [Tl] we introduce two exact, non-perturbative 
univariate distributions which produce non-Gaussian ensembles with fixed variance and 
skewness. In section |TJ we derive the higher moments for both distributions. Section [TV| 
shows how we can employ the non-Gaussian 1-D distributions to generate isotropic en- 
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sembles of non-Gaussian maps with given angular power spectra and bispectra. We show 
examples of such distributions in section |V| and discuss the results and applications of this 
work section |VT[ 

II. EXACT NON-GAUSSIAN 1-D DISTRIBUTIONS 

In this section we introduce two classes of distributions which can be employed to generate 
random numbers with exactly specified second and third moments. The first was introduced 
as a non-perturbative, non-Gaussian likelihood in the Bayesian analysis of the CAT/VSA 
data 0. It was originally introduced, in a theoretical context, in the study of off-the-ground 
state perturbations in the inflaton field [§]. The second class of distributions is an ad hoc 
solution obtained by requiring the simplest superposition of Gaussian PDFs that results in 
a distribution with a given second and third moments. 

Both PDF functions can be used to generate distributions with fixed power spectra 
and bispectra (and optionally higher order spectra). The first is best suited for mild non- 
Gaussianity; the latter for strong departures from Gaussianity. 



A. The non-perturbative harmonic likelihood 

We now summarise the method of Rocha et al. [0J§. Let x represent a general ran- 
dom variable. We build its distribution from the space of wave-functions which are energy 
eigenmodes of a linear harmonic oscillator (see e.g [10]). We have that 

i>i x ) = ^2 a n^n(x), (1) 
where n labels the energy spectrum E n = hu(n + 1/2), and 

^(aO = tf»#«(;^)e~^, (2) 
with normalisation fixing C n as 
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The only constraint upon the amplitudes a n is 

EKI 2 = i, ( 4 ) 

which can be eliminated explicitly by imposing the condition 
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Oo is the variance associated with the (Gaussian) probability distribution for the ground 
state iV'ol 2 - We define Hermite polynomials H n {x) as 

Hn (x) = (-lTe* 2 £;e-* 2 , (6) 

with normalisation 

/oo „ 
e~ x H n (x)H m (x)dx = 2 n ^n\5 nm . (7) 
-oo 

The most general probability density is thus: 

_ * 2 

The ground state (or "zero-point" ) fluctuations are Gaussian, but any admixture with other 
states will be reflected in a non-Gaussian distribution function. 

It can be shown ||[7| that the a n reduce to the cumulants K n (up to a multiplicative 
constant) for mild, "perturbative" , non-Gaussianity. This is achieved by reducing the prob- 
ability density (|8]) to an asymptotic expansion around the Gaussian distribution function 
(the so-called Edgeworth expansion). Such asymptotic expansions on the space of orthonor- 
mal Hermite polynomials suffer from the fact that truncations at a finite order of cumulants 
lead to pseudo-distributions which are not positive definite. This is not the case here and 
the advantage of using the a n over cumulants is that setting all but a finite number of 
them to zero still leads to proper distributions. In some sense the a n are non-perturbative 
generalizations of cumulants. 



The above probability density may be easily applied to generate a centered distribution 
with a fixed variance and skewness. Let us start with the PDF 
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a + -7=^1 -tt" + ~^H 2 —f=- + -r^H 3 



, (9) 



in which all the a n are real and we have set all a n = forn > 3. We then calculate moments 
around the origin defined as 

A(n = ( x n ) = r x n P{x)dx. (10) 

<J — oo 

For a normalised density we have /x = a 2 , + a l + a 2 + a 3 = 1 an d the first three moments 
given by 

//i = (2o- 2 ) 1/2 (^2q; 1 q;2 + v^aoai + VQa 2 a 3 ^ , 

A*2 = o" 2 («o + 3a 2 + 5a 2 . + 7a 2 + 2^0:0^2 + 2v / 6o; 1 a3^ , 

(3 9^/3 \ 

-j=a ai + 6aia 2 + ~^-a 2 a 3 + VSa a 3 \ . (11) 

The aim is to have a centered distribution with jii = 0. This can be achieved by setting 
«i = and a 2 = in which case we are left with the system 

1 = al + a 2 , 
/i 2 = (T 2 {al + 7a 2 ) , 

^3 = (2a 2 ) 3 / 2 v / 3a «3- (12) 



Setting a = V 1 - a 2 , we then solve for a and 0:3 for particular values of the required second 
and third moments ji 2 and /i 3 



2V6(1 -ai) 1 / 2 ^' 

a problem we deal with numerically. Hence we arrive at a 1-D centered distribution with 
fixed variance and skewness. 



X X 

FIG. 1. The left panel shows the distribution P{x) for various values for the relative skew- 
ness parameter s = fis/fj^ (0.0 < s < 0.7) . The right panel shows the distribution T(x) for 
0.0 < s < 10.0. 

By restricting ourselves to only two parameters we have in fact constrained the space 
of possible distribution functions. This is reflected in the fact that we cannot generate 
distributions with any given variance and skewness. By studying the function for the relative 
skewness s 

_/i3_ 2y/E(l - a 2 3 )a 3 . . 

(l + 6a§)»/* U4J 

it is easy to see that the maximal relative skewness is 

s = ±0.739 at ct 3 = ±0.278. (15) 

In general our method can generate higher values of s (since it can generate any distribu- 
tion), but for that purpose one needs more parameters a n . The generalization of the above 
construction for more parameters is trivial if somewhat tedious. 



6 



B. A tailor made distribution function 



One way to overcome the limitations of the distribution described above without increas- 
ing the number of parameters is to build a distribution by 'summing' Gaussians. In fact it 
is obvious that the distribution P(x) can be well approximated by a series of uncentered 
Gaussian functions of different heights and widths. Consider, for example, the distribution 
obtained by the superposition of three Gaussian distributions 

T(x) = — L ( e - 2 + c -C-AA) a + (3 e-^ 2 ) . (16) 
(2 + AOv 71 " 

In this case the first three moments of the distribution are given by 



m = o, 



L- T (2 + ft + 2fl ) /3 1 2 (l + A))) 



^"2(2 + A, 

" S - 2 + (3 ■ (17) 
In contrast to the previous distribution the value of s is now effectively unconstrained. In 
practice we find that in this case the range of s is limited only by numerical effects in the 
random number generation process. 

Fig. [l] compares the distributions P(x) and T(x) for various values of the relative skewness 
s. The former may be considered somewhat 'unphysical' since the minima in the functions 
will tend to produce 'holes' in the distribution of variates for each particular set of moments 
(this is not surprising since the PDF originates from the wavefunction of an oscillator). We 
should point out though that the present exercise is aimed at producing maps for the sole 
purpose of testing statistical tools and not to reproduce theoretically motivated models of 
the cosmic microwave backgrounds. The adaptation of the method to more suitably physical 
non-Gaussian distribution functions such as the PDF T(x) above is trivial as we have shown. 

A word of caution is in order, concerning the numerical aspects of generating univariate 
random numbers with a given non- Gaussian PDF. The Jacobian required to generate the 
random variables using the transformation method is not easily computable in the case of 



the PDFs P(x) and T(x). We advocate, instead, the use of the well known rejection method 
|IIJ to generate the required one dimensional random variables. 



III. LEAKAGE INTO THE HIGHER MOMENTS 



Both distributions have general solutions for all higher order moments fi n . The fourth and 
sixth order moments are of particular interest since they generate the sample variance and 
hence the cosmic variance in the observed power spectrum and bispectrum of the realisations. 
For n even we find the following relatively simple expressions for the moments 



fJ-n 



2-n 
«0 F 



n 



+ ai{ 3r 



3 + n 
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5 + n 



4 /7 + n 
3 V 2 



(18) 



for the distribution P(x) and 

r(^) 



/in 



(2 + /? )xAF 



e f 1±^, I ft) + e-^M (1±±, + 1 



(19) 



for the distribution T(x) where M(a, b, c) are the confluent hypergeometric functions |T9| 



Unfortunately the non-linear dependence of the parameters a, a n and (3 n on the required 
ensemble moments complicates the analytical form for the sample variances u 2 ^^) and 
a 2 (fi^). In practice therefore the sample variances would be calculated numerically via the 
above expressions and making use of the solutions to the systems ( |12| ) and ( |TTD or via 
Monte-Carlo simulations. 



IV. GENERATING NON-GAUSSIAN MAPS WITH FIXED ANGULAR POWER 

SPECTRUM AND BISPECTRUM 

We now propose to apply this method to the generation of non-Gaussian maps with fixed 
angular power spectrum Ce and bispectrum Bf . The broader context is the generation of 
non-Gaussian signals to be be used in simulations of upcoming satellite experiments. 

The distributions of the previous section have the disadvantage that they can only gen- 
erate non-Gaussian 1-D distributions. As a result, when a likelihood constructed using P(x) 
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is applied to CMB maps as in Rocha et al., it probes a combination of non-Gaussianity 
and anisotropy. Consider, for instance, the generation of a full-sky map with temperature 
fluctuations ^(n), by means of its spherical harmonic components: 

AT 

— (n) = ^2a em Y em (n). (20) 

If we set «3 7^ in the above distribution and use it to generate a set of 3l(aeo) we will 
have ((JRctio) 3 ) ^ 0, with ((Jftae m ) 3 ) = for m / 0. Isotropy, on the other hand, imposes 
"selection rules" upon correlators, in this case 



a hm 3 ) 



f tl h £3 ^ 



BtiMsi (21) 



y mi m 2 m 3 J 

where the quantity (. . .) is the Wigner 3J symbol, and the coefficients Be^h are the bis- 
pectrum (abbreviated to Be for the case £\ = £2 = £3)- Hence the distribution we have just 
generated is not only non-Gaussian but also automatically anisotropic since all third order 
correlators except for (a^o) are zero. 

If we are to impose isotropy, we must necessarily have correlated ae m , a feature not 
allowed by the method of Rocha et al. In the present context a correlated set of ag m is 
equivalent to a series of coupled harmonic oscillators. The obvious way to achieve the 
necessary correlations would be to use the Hilbert space of coupled harmonic oscillators to 
set up the most general multi-variate distribution but this proves to be impractical. Here 
we introduce a simple but crucial modification which restores isotropy. 

We propose to set up an isotropic ensemble in two steps. First we generate an anisotropic 
ensemble by drawing all ae m from a Gaussian distribution with variance spectrum Ce, except 
for the m = mode. The latter is given one of the PDFs described in the previous section, 
with variance Ce and skewness 

/ V 1 

000 



B e (22) 



y0 y 



The appropriate PDF is enforced using the rejection method, as described above. For each 
set of Ce and Be we solve the systems given by Eqs. fljjD or ([FTP numerically to obtained 
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the required values for the parameters a, ao and 0:3 or (3 and Pi . We label the resulting 
ensemble a, for anisotropic. 

We then apply a random rotation to each realization in the ensemble a, with a uniform 
distribution of Euler angles. In this way we arrive at an isotropic ensemble of temperature 
maps, labelled by i, with spherical harmonic coefficients 



*4 = £^L^)<4 



(23) 



where is the rotation matrix, and Q denote the 3 Euler angles. The 2-point correlators 
are 



(24) 



hence we have CI = Cf. 



The 3-point correlators are zero for any correlator involving different I. However we now 
have 

= (KAmLornKsoi^n >< s * 
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(25) 



mi m 2 m 3 

(note that with one index m set to zero the rotation matrices reduce to spherical harmon- 
ics, one of the Euler angles becoming irrelevant). Again we have recovered the original 
bispectrum with B\ = B%. 

Hence the procedure we have defined produces the desired isotropic ensemble with a fixed 
angular power spectrum and bispectrum. The random rotations produce the necessary 
correlations between the a^ m coefficients to ensure isotropy. However, by means of this 
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procedure, we are unable to generate maps with non vanishing inter-£ bispectrum coefficients 
(as studied in 

Note that we can also draw all ai m from the 1-D PDFs introduced above, and then 
subject them to a random rotation. The argument presented here still goes through, and 
an isotropic ensemble is still obtained. However, as we shall see later, the cosmic variance 
in the estimators for the Ci and Bi will be larger in this case. 



V. SIMULATED NON-GAUSSIAN MAPS 

We now present CMB temperature maps using the prescriptions detailed in the previous 
sections. We generate full sky maps pixelised using the HEALPIX 0] package at pixelisation 
levels 128 and 512, equivalent to 196, 608 and 3, 145, 728 pixels respectively. We generate 
maps which include a cosmological signal in the form of a standard ACDM power spec- 
trum computed using CMBFAST jD]] and finite beam sizes. Instrument specific noise and 
physically motivated profiles for the bispectrum will be treated in further work ||15| . 

For simplicity we shall assume a bispectrum which is fixed relative to the power spectrum 
(i.e. Bi = sC^ 2 ) with a ratio which is the same for all scales (i.e. with S£ = s constant). 
This is a rather artificial asumption since in practice the host of non-linear, secondary and 
foreground contributions to the CMB (T^| and possible primordial non-Gaussianity will arise 
at different angular scales. Theoretically motivated bispectra and the addition of instrument 



specific noise will be discussed in further work [15|, and are straightforward to obtain with our 



method at no extra computational cost. We use both the PDFs discussed above to generate 
the 1-D distributions required for the production of the initial anisotropic ensembles. 

In Fig. Q we show a Gaussian realisation of a standard ACDM model at COBE-DMR 
PDfl resolution. The equivalent non-Gaussian realisations with the same power spectrum are 
shown in Figs. ^| and |] which are generated using the PDF T(x) and values for the parameter 
s of 0.5 and 2.0 respectively. Fig. |5| shows a non-Gaussian realisation (s = 2.0) with the 
same power spectrum and MAP [^TJ resolution. 



11 




FIG. 2. HEALPIX map of a noiseless, Gaussian CMB realisation at COBE beam resolution 
(normalised temperature units). The map was generated using a standard ACDM power spectrum 
(J2a = 0.7, Qcdm = 0.25, Q b = 0.05, h = 0.7 and n s = 1.0). 




FIG. 3. A non-Gaussian realisation at the same resolution generated with the same ACDM power 
spectrum and a 'white' bispectrum with s = 0.5. The PDF T(x) was used as the non-Gaussian 
generator. 
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We first make use of our maps to examine possible effects of non-Gaussianity upon 
estimation of the power spectrum. For this purpose we compare the observed power 

1 



a 



obs x l„^2 



I E Kl, (26) 



2i+ , 

m=—t 



in the non-Gaussian ensembles with that of the usual Gaussian ensembles. Different PDFs 
will produce ensembles with different cosmic variances (for the power) and in particular the 
cosmic variance will be different from that of Gaussian ensembles. Assuming the a m modes 
to be Gaussian distributed the power spectrum will have a xli+i distribution which yields 
a particularly simple sample variance 

ACt) - (27) 



for zero noise and infinitely thin beam (see e.g. |12|,|T8[). In the non-Gaussian case the 



variance will assume a more complicated form due to a non-zero contribution from the fourth 
order cumulant (or connected moments in analogy to the connected Green's functions of field 
theory) : 

cr 2 (n2) = K 4 + 2/i 2 . (28) 

In many motivated examples non-Gaussian processes display higher cosmic variance in the 
power; an example is texture models as studied in |E| . 

In Fig. |6] we show the average power spectrum for an ensemble of 10 maps for mild and 
extreme values for the relative skewness parameter s. The maps do not include any noise 
contribution or beam effects. The right panel shows the average power for maps where all 
a£ m modes are non-Gaussian in the original anisotropic, uncorrelated ensemble. The left 
panel instead is for the case where only the a^ of the original modes are non-Gaussian. 
The dotted line shows the extent of the 1-sigma cosmic variance for an equivalent Gaussian 
ensemble. 

We see that the sample variance of the power for the non-Gaussian variables is compa- 
rable to that of its Gaussian counterpart for low s and grows for larger value of the relative 
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FIG. 6. The estimated power spectrum for ensembles of 10 pure signal maps. The solid line is 
the power spectrum of a standard ACDM model (Oa = 0.7, tie dm = 0.25, = 0.05, h = 0.7 
and n s = 1.0). The dotted lines show the extent of cosmic variance of the averaged power for 
the respective Gaussian ensemble. The short dashed line is for an ensemble generated using the 
distribution P(x) with s = 0.5. The short dashed dotted line is for s = 0.5 but generated using 
T{x) and the long dashed dotted line is for T(x) with s = 5.0. The ensembles in the left panel 
were generated using only the a^Q modes as the non-Gaussian seeds in the original anisotropic, 
uncorrelated ensembles whereas in the right panel all the modes were originally non-Gaussian. 
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skewness. This dependence is expected since we know that k 4 (ct, «o,«3) = ^(s, Cg). Note 
that the isotropising rotation does not affect ensemble averages, so that the sample variances 
of estimators in the isotropic ensemble trace the variances of the 1-D distributions generat- 
ing the non-Gaussianity. In general the contribution from the connected moments can be 
derived analytically since they can be expanded in a series of moments |Tj|. In this case 
though the solution for u{^2) in terms of Cg is very contrived and has no simple expression 
as the one above for the Gaussian case. Numerically estimation on a case by case basis will 
give a more feasible determination of the sample variance of the estimators being used. 

We now turn to the issue of the detectability of non-Gaussian bispectra making use of 
our maps. In Figs. ^ and [| we show the distribution of the value for the observed bispectrum 



( \ 

a m 1 a mz a m s i (^9) 

i mi Tn 2 m 3 



as measured from 20000 realisation ensembles. Our Wigner coefficients are calculated using 
the recursive relations of Schulten et al. [14j] and are accurate to £ > 3000. 



We show the histograms for multipoles £ = 2,4,6 and 8 and for values of s — 2.0 and 
0.5. We used the distribution T(x) to assign non-Gaussian values to all the ag m modes. 
It is interesting to note that the sample variance of the estimator B t hs is reduced in the 
non-Gaussian case with respect to that of the Gaussian ensemble. This shows that the 
contribution from the connected moments to the non-Gaussian part of /xg is negative and in 
particular 

|k 6 | > 15K 4 /i 2 + 9«3, (30) 

confirming that the leakage into higher orders is not negligible for both PDFs. 

The implications of this result are interesting. It seems that, for the type of non- 
Gaussianity which we have generated, when searching for non-Gaussianity and armed with 
a single realisation, we could only rely on the rejection of the non-Gaussian hypothesis by 
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FIG. 7. Histograms for the measured value of the estimator B° for 20000 realisations. The 
solid line corresponds to the non-Gaussian ensemble generated using the distribution T{x) with a 
value of s = 0.5 and the dashed line corresponds to the Gaussian ensemble. 

observing modes whose bispectra are (counter-intuitively) too large for them to be non- 
Gaussian. Conversely if we were to measure bispectra which are consistently too close to 
zero for enough I we could reject Gaussianity in favour of non-Gaussianity. Even though this 
remark may sound counter-intuitive it is the basis of the result in || (in which Gaussianity 
was rejected on the basis of a low x 2 ). 

For completeness we have displayed in Table | all the results concerning C^ bs and B^ bs 
in both the Gaussian and non-Gaussian ensembles as discussed in the previous paragraphs. 
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2 


892.729 


890.716 


892.590 


-4938.244 


-325.169 


-6376.174 


4 


259.111 


258.795 


257.526 


426.451 


-97.583 


559.303 


6 


120.751 


120.684 


121.219 


-160.554 


-4.521 


-123.481 


8 


70.419 


70.401 


70.388 


13.399 


-11.157 


42.092 


10 


46.657 


46.656 


46.687 


-27.442 


3.777 


-18.385 


12 


33.704 


33.793 


33.783 


5.682 


-8.440 


9.484 


14 


25.713 


25.674 


25.763 


-6.365 


-5.433 


-5.449 


16 


20.419 


20.480 


20.383 


1.775 


0.546 


3.389 


18 


16.747 


16.737 


16.730 


-4.908 


-0.757 


-2.245 


20 


14.084 


14.087 


14.077 


2.081 


0.431 


1.563 



TABLE I. The observed statistics for the first 10 even multipoles for both a Gaussian and 
non-Gaussian noisless ensemble. 

VI. DISCUSSION 

In this paper we proposed a method with which to generate non-Gaussian maps with 
fixed power spectrum and bispectrum. Our strategy is as follows. We generate these maps in 
harmonic space, and give each of the modes, independently, a 1-D non-Gaussian PDF. The 
resulting maps are anisotropic, but a random rotation then restores the ensemble isotropy. 
We proved that the power spectrum and bispectrum of the isotropic ensemble is simply 
related to the variance and skewness of the 1-D PDFs employed. We then proposed two 
different PDFs for which the variance and skewness may be fixed independently. One is 
based upon the Hilbert space of an harmonic oscillator; the other upon a superposition 
of Gaussian functions. In both cases we worked out the algebra relating the parameters 
controlling the distributions and the variance and skewness. We then drew our random 
numbers numerically using a rejection method. 
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Even though we only considered the generation of maps with a fixed power spectrum 
and bispectrum we stress that the extension of our method for higher order moments is 
straightforward, even though it has its computational costs. The various parameters of the 
PDFs considered may be used to fix the higher order cumulants of the mother 1-D PDF 
applied to single modes. Once again a random rotation restores isotropy, and the isotropic 
higher-order correlators may be simply related to the cumulants of the original PDFs. 

We close with an evaluation of the efficiency of our method when meeting the ever 
improving resolutions of upcoming experiments. As we have shown it is possible to generate 
high resolution non-Gaussian maps with this method. However they are computationally 
more expensive to produce than Gaussian maps, specifically due to the random rotation to 
which they must be subject. Currently the rotation is carried out in £-space which requires 



the computation of the rotation matrix elements |23| . This avoids sampling problems which 
would arise if it were carried out in pixel but requires increasingly long computation times for 
increasing I. On the other hand the real bottleneck is likely to be at the analysis stage, rather 
than in the simulation of maps. Sums involving Wigner coefficients become computationally 
intensive as i is increased. 
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